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ABSTRACT 


This paper studies alternate ways to model density discontinuity in split-step Fourier 
parabolic equation models. The Monterey-Miami Parabolic Equation model is used to 
implement an alternative to the effective index term in the smoothing function and a split- 
step Fourier/Finite Difference hybrid technique. The model is shown to converge to a 
stable solution that is slightly lower than the benchmark solution. A range step size of 
approximately one wavelength is shown to provide the closest approximation to the 
benchmark solution. Acceptable solutions are obtained with large depth grid sizes for the 
alternate smoothing function. Smaller depth grid sizes are necessary for accurate 
solutions when using the hybrid implementation technique. The effect of reference sound 
speed is shown to minimize the phase error present when the models are used in the 
presence of a strong density discontinuity. 
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I. 


INTRODUCTION 


A. BACKGROUND 

Parabolic equation (PE) methods have been used in the field of physics for many 
years to reduce the two-way, elliptic wave equation to a simpler, one-way equation. 
Although originally utilized in the field of plasma physics, Hardin and Tappert [1] were 
the first to introduce the PE method for solving underwater acoustic propagation in the 
early 1970s. The PE approximation to the wave equation is first order differential in 
range (r), and therefore was able to be numerically solved by using a one-way marching 
algorithm, such as the split-step Eourier (SSE) algorithm. The SSE technique works by 
first sampling the acoustic field on a depth grid, i.e., with uniform sampling Az, at a 
single range r. The depth grid then undergoes two fast Eourier transforms (EETs) and 
two vector multiplications in order to obtain the solution at the next range step r-l-Ar. 
The process is repeated at each range step until the entire acoustic field is modeled [2]. 
The major advantages of the SSE technique are its computational efficiency, stability, and 
simplicity. 

Alternately, finite-difference (ED) techniques may be used to solve PEs. 
Numerical implementation of the ED algorithm is similar to SSE. The acoustic field is 
first specified at all depths at an initial range. A finite-difference approximation based on 
a Taylor series expansion of the derivative terms is then directly applied to the density- 
dependent, transverse differential operator. The solution is obtained by solving a 
tridiagonal system of equations with dimensions equal to the number of depth samples at 
each range step [3]. By directly applying the finite-difference approximation to each 
element in the depth mesh, the ED technique avoids having to use a smoothing function 
and can be used in a variety of ocean environments but at the cost of increased 
computational complexity. 

A third method of solving PEs is the finite-element (EE) method. This method 
discretizes the ocean environment into multiple nodes using Galerkin’s method [4]. Then 
by applying continuity of pressure and the normal component of particle velocity to the 
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boundary conditions of each node, a tridiagonal system of equations is established that 
can be solved through well-defined methods. By discretizing the environment, the FE 
method has the ability to increase or decrease the number of nodes within a local depth 
range, which leads to more detailed or faster calculations than possible by using FD 
techniques [5]. However, due to the complexity of FE schemes and their associated 
stability issues, FD and SSF techniques are the preferred method in the modeling 
community [6]. 

The 1981 PE Workshop [7] proved Harding and Tappert’s SSF technique was 
able to efficiently and accurately model the acoustic field in long-range, narrow-angle 
problems with negligible bottom interactions. More advanced SSF methods based on 
higher-order operator splittings of the so-called “square root operator,” such as 
introduced by Thomson and Chapman [8], also greatly improved the accuracy of the 
approach for wide-angle solutions. 

In areas with strong bottom interactions, such as in shallow water, the SSF 
technique has been shown to be less accurate. This is due to the fact that the propagator 
term in the SSF algorithm contains density derivatives in the effective index of refraction 
term, which require the use of a smoothing function to approximate the density profile in 
the vicinity of the ocean bottom [9]. This accuracy of the smoothing function’s 
approximation is considered to be the major weakness to SSF based models. 

In 1997, Yevick and Thomson [10] proposed a hybrid SSF/FD method. Their 
initial motivation was to introduce a small set of equations (generally five) that treat the 
density discontinuity directly through a finite difference approximation rather than 
relying on the mixing function previously employed, thereby improving the accuracy of 
the solution. The remainder of the calculations and marching algorithm are performed via 
SSF techniques. This method capitalizes on the accuracy of the FD technique while 
utilizing the efficiency of a SSF marching algorithm. This method has shown accurate 
results and is efficient enough for use in two and three dimensions. 

In addition, they were able to show some improved phase accuracy (through less 
sensitivity to the reference sound speed) by incorporating additional higher order terms in 
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the operator splitting of the square root operator. These additional terms were also shown 
to lend themselves to finite difference solutions. In this case, however, the application of 
the set of finite difference equations extended over the entire depth computational 
domain. 

B. PROBLEM STATEMENT 

The Monterey-Miami Parabolic Equation (MMPE) [11] method utilizes a SSE 
technique to model the propagation of the acoustic field underwater. Versions of MMPE 
have been developed for both 2-D and 3-D, single frequency and broadband calculations, 
and including various types of environmental fluctuations, such as internal waves, surface 
roughness, 3-D bathymetric variations, turbulence, and others. The MMPE model utili z es 
the wide-angle PE (WAPE) approximation of Thomson and Chapman [8], and deals with 
the density discontinuity by using a hyperbolic tangent mixing function approximation to 
the Heaviside step function at the bottom interface. 

This model has been shown to be accurate in deep-ocean problems [12]. 
However, MMPE’s treatment of the bottom condition can become problematic in shallow 
water environments where bottom interactions play a significant role. In this thesis, we 
will consider solutions to the acoustic propagation in a simple, Pekeris waveguide, where 
the environment is defined as having a lossless water column with constant sound speed 
of 1500 m/s and a bottom depth of 300 m, below which is a homogeneous half-space of 
sound speed 1700 m/s, density 1.5 g/cm3, and attenuation 0.5 dB/A. The source is 
modeled as an omnidirectional point source at depth 180 m transmitting a single- 
frequency, continuous wave (CW) signal at 100 Hz. Eor comparison purposes, a 
transmission loss trace will be extracted at 100 m depth out to a maximum range of 
20 km. 

When compared to a coupled-mode model, known as Couple (and considered 
capable of generating benchmark quality solutions) and a finite-element parabolic 
equation (EEPE) model, the MMPE model is quite accurate at shorter ranges, as can be 
seen in Eigure 1. However, at about 5 km a phase error begins to appear and accumulate 
as the range is increased. The objective of this work is to study multiple implementation 
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options for modeling the density diseontinuity in MMPE in order to improve the model’s 
long-range aeeuraey. Insights into the eause of this phase error aeeumulation with range 
will also be explored. 


tawh-iho()•*» iQDODIODlim lODODft 16D0D* tOD 





Figure 1. Plot of TL (dB) vs. Range for Couple, FEPE, and MMPE. 
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II. THEORY 


A. PARABOLIC EQUATION DERIVATION 

We begin with the Helmholtz equation for acoustic pressure, p(^r,z), in 
cylindrical coordinates for a time-dependent harmonic field, 


lA 

r dr 


dp 

dr 




1 dp 
p{z) dz 


+ kln^{r,z)p = Q, 


( 1 ) 


where kr .=— is the reference wavenumber, n(r,z) = —-— is the acoustic index of 
Co c(r, z) 

refraction, Cq is the reference sound speed, c(r,z) is the acoustic sound speed, and p(z) 

is the density (which is assumed to vary only with depth at any given range). It is within 

c(r, z) and p(z) that all features of the environment are represented. 

If we now introduce the so-called PE field function according to 


^(r,z) = -^p(r,z) 




( 2 ) 


and assume that the solution is dominated by the forward propagating field, it can be 
shown [9] that the 2-D Helmholtz equation can be reduced to a parabolic equation of the 
form 


dr 


a,(e-i)'P , 


(3) 


where the square-root operator, Q , is defined as 


By defining the operators 


Q = 


2 l d 


11 

pdz) 


(4) 
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(5) 


s = n^ — \ 

^ kl dz^ 

_ I dp d 
^ kip dz dz 

the square-root operator may be written as 

Q = ^\ + £ + p + y . (6) 

Solutions to the parabolic equation may then be obtained from a one-way 
marching algorithm solution in the general form 

'F(r-rAr,z) = exp[jloAr(g-l)]'P(r,z) . (7) 

In practice, it is not possible to directly apply the square-root operator in the marching 
algorithm defined above. Some form of approximation is required. 

Alternatively, we could define the field function as 

^ = yfp^, ( 8 ) 

which then satisfies a PE of similar form 


In this case, the square-root operator is defined by 

where 

£ = n^ -1, and =n^ ^ 

2kl 


I gy 3^1 dp\ 

p dz 2yp dz ^ 


(9) 


( 10 ) 


( 11 ) 
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The original approaches in MMPE were based on this latter definition, which requires a 
description of p{z) in terms of a mixing function to smooth out the discontinuity at the 
bottom interface so that the derivatives may be computed. The solution to this PE is also 
obtained from a marching algorithm of the form given in Equation (7). 

B. APPROXIMATION OF THE SQUARE-ROOT OPERATOR 

The Thomson-Chapman WAPE operator splitting [8] produces an approximation 
to the square-root operator in Equation (6), defined by [10] 

= Vl+£■+^ • (12) 

By using this definition, the propagator in the marching algorithm becomes 


exp [z^oAr (21 -1)] = exp id ^ ^/lT^ + yjl + ju -f- ^Jl + / - 3 j 

ij exp 


« exp 

1 1 

Si 

-1- 

1 

1 _ 1 

exp 

zA(7i + /^-1) 

= exp 

zA(7i + /-1) 

exp 



1 + s-l 


(13) 


where S = k^Ar, and the approximation in the line is due to the neglect of any non¬ 
commutation of the operators £, p, and y. As noted by Yevick and Thomson [10], the 
order in which the operations are implemented is important, and the term with y should 
be applied at the end of the range step. Note that the operator currently used in MMPE 
corresponds to 

2i = Vl + ^ + \l\ + ju — \ , (14) 


and there is no y term of concern. 


The Thomson-Chapman approximation neglects the lowest order cross terms 
between the operators. The lowest order correction is [10] 


22 

= Q'i-\'i£{M+r)+{M+r)£+My+m'\ 

o 


(15) 
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For completeness, note that the corresponding lowest order correction to the existing 
MMPE operator would be 


e;=e;-- 



8 


[^/y + /y£^] 


(16) 


In PE models based on the SSE algorithm, the operators containing S (or s) and 
// are treated as simply scalar operators in the physical space (z) and wavenumber space 

respectively. The vector multiplication is then a trivial task, which leads to the 


great efficiency of the SSE approach. However, this does require that operators 
containing s and // are separated in the approximation of the square-root operator. 


If we assume that the y terms and cross terms can be added later as corrections, 
then an intermediate step solution can be defined as 




r 

f 1 7 2 ^ 


T'(r + Ar,z) = EEr ‘ 

expj 



>xFFT |exp[/A(n-l)]'P(r,z)| 


Note that if we replace nhy h, then this is the complete range step utilized in the existing 
MMPE code. 


C. DENSITY SMOOTHING APPROACH 

In order to utilize the simpler approach in MMPE, based upon Equations (9), (10), 
(11), and (14), we must first define the mixing function for the density contrast. In this 
case, the density profile is defined in terms of a generalized function representing the 
discontinuity at the water-bottom interface by 

P(z) = p„+(Ph-p„)H(z-z^), (18) 


where 



1-i-tanh 

fO’ 

'' ’ 2 


[2l)_ 


(19) 
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is the hyperbolic tangent mixing function approximation to the Heaviside step function. 
The parameter g = z — Zi,, where Z^ is the bottom depth, and L is the mixing length of this 
smoothing function. 

Now, since = n^ -^ 

Ikl 


1 ay 

p dz 


^ 1 dp^ 
p dz 


, we need to define 


dz dz 


AL 


\2Lj 


( 20 ) 


and 


aV . . d^H 

= [Ph-P«y 


{Pb-P«) 


dz 


dz' 


AU 


tanh 


\2Lj 


r 


sech^ 


\2Lj 


( 21 ) 


Previous analysis by others [7], [9], [11], [12] has shown that good solutions result by 

2 ;i 

choosing a mixing length of L = — ~ . Although this is an additional free parameter 

K ^ 

that could be evaluated for accuracy, we shall not explore this further in this thesis and 
will keep this mixing length defined throughout for those algorithms that utilize the 
density smoothing approach. 


D. YEVICK AND THOMSON HYBRID IMPLEMENTATION TECHNIQUE 

The hybrid approach by Yevick and Thomson [10] shows that the new formalism 
defining Q[ requires evaluation of the additional operation 


^'(r +Ar,z) =exp + y-1^ 'P(r + Ar,z) 


( 22 ) 


Note that this completes the range step if we only invoke the QI approximation, and we 

may replace T*' ^ ^. The Q 2 approximation requires an additional step and is 
discussed later. 

Yevick and Thomson [10] proposed treatment of this operation through a finite 
difference approach, whereby the Fade approximation 
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(23) 


exp 


l + -{l + iS)/ 




generates the finite difference equations 


1 + ^( 1 -^'^)^ 


'F(r + Ar,z) = 


\+^{\+i5)r 


'F(r + Ar,z) . 


(24) 


dp 

In this case, direct application of y is problematic since y depends on —, and 

dz 

we do not want to invoke the smoothing function for p . Yevick and Thomson [10] show 
that an alternative is to evaluate a finite difference approximation to the original 
transverse operator 


kl dz 


p dz 


(25) 


Then y - p' - p , where p is well-defined as an operator on 'P . 


Specifically, using a 3-point centered difference approximation defined by 

( \ 


f ~ — 

Jo , 
h 


/i-/i 

V 2 2 J 


+ 




(26) 


then 



where the indices 0, ± —, ±1 apply to the depth index in question and neighboring 
depths above and below. 
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Note that Yevick and Thomson assumed a eontinuous funetion for p , sueh that 


A =^(a+a) > P 1 =\{p^ + P-x) • 

2 ^ 2 ^ 

However, we shall treat the density as defined by 


P{z) = P^+{Pb-P.)^{z-Zb) = \ 


Pw , z < Zf,, 

|(a+a) ,z = zp 


Pb 


,Z>Zu. 


(28) 


(29) 


Then, sinee 


kl dz kl dz V dz 


k^Az 




dz 


2 


az 


2 


(30) 


1 


(koAz)' 




(kyAz)' 


-['T,- 2 'Ty+'T_J 


we obtain 


f¥ = {p'-p)W 




r ^ 



Po _ 2 

'I^i- 


Pi 



L 2 J 

V 



Pi 


2 


_P^ 

P 1 



1 1 

r ^ 


2 

o 

+ 

Po 2 

^ 1 

1 






1 2 J 



• (31) 


This agrees with Yeviek and Thomson when their notation is employed, defined as 


A 

Pi 


Po 


(a+Po) 


= P+ 


Po _ A) 

\{p, + P-i) 


Po I Po 

Pi P i 

2 2 


= 2Po 


- + 


(Pl+Po) (Po+P-l) 


= P++P- 


(32) 
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For our discontinuous density, only integer index values of p have any meaning. So we 
shall define the notation p^ = p^= p{z> ■, and p j = p_= p{z<z^). 


Then, 


f¥ = 


1 


(k^Az) 




( 

\ 


f \ 

” 

Pq _ 

'i'l- 

Pq I Pq 

-2 

'i'o + 

Pq _ 

'i'-i 

_V/^+ y 


kP^ P- 

/ 





(33) 


We can then use a tridiagonal matrix solution to the expression 

41^ = 5^, 


(34) 


where A has elements defined by 


1 


1 + 


and B has elements defined by 


i+-{i-is)r 


. Note that A and B are just identity matrices except in the vicinity of 
Z~ Zj, , where y is non-zero. 


If we are only interested in the correction from the definition of , the solution 

is complete. However, Yevick and Thomson also go on to describe how this can be 
further refined through application of another finite difference approximation to the 

additional term in ^,i.e., +y^ + {^p + y^£ + py+ yp^. Since p + y^p’, 


we may write this as 


1 


r 1 \-—i5{£p' + p'£ + py + yp) 

expj-—[£•//'+ //'^+ /// + ///] U—- 

' 1-1-— iS [£p'+ p'£ + py + yp') 

16 


(35) 


Note that for y - 0, this correction simplifies to 


1 


f ,1 + 

exp ^-—^ - , 

^ l+—iS[£p + p£) 


(36) 
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which can be applied as a eorreetion to the existing MMPE approaeh based on Q 2 if we 
replaee s^s. 

Yeviek and Thomson drop the /// and yju terms, arguing that these are small, 
higher-order eorreetions, and instead utilize the eorreetion [10] 


exp 


y 


+ /u's) 




l + ^idi^s/u' + fj!s) 


(37) 


This produces the additional eorreetion for 
'T (r -I- Ar, z) = exp 




'Fy-I-Ar, z) 


(38) 


where we again invoke the matrix approaeh 

(39) 


and elements of A are defined by l+—iS(^£ju' + ju'£) while the elements of B are 
defined by l- — iS(^£ju' + ju'£). 

The evaluation of f/zT follows Eq. (27), 


£jU'^ 


^oPo 


(koAz) 



r 3 


1 

1 1 

1 


— -1- — 

vp -1- vp 

^ 0 ^ ^ -1 

p\. 

: Pi p 1 : 

pj. 

2 

(2 2 ) 

2 


(40) 
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However, for jj’s^ we must consider 



"i 

A 5 

■j_r 

kl dz 

p dz 

^0 Az dz 



Po 


(^qAz) 


V 2 2 2 2j 


Px P 1 


(41) 


Po 


(^o^z) 



f ^ 


1 ... 

1 1 

1 


—+ — 

^o'Po + ^-i^'i'-i 

Pl 

Pi P 1 

Pi 

2 

1 2 1 ) 

2 


Note that if we are applying only this correction to MMPE, then the correction 
would result from the application of Equation (36) to A'E(r +Ar,z) = B'E'(r + Ar,z). 
Here, 


efP¥ « 






(42) 


and 


juP¥ « 





(43) 


As before, we would need to replace if we are utilizing the density smoothing 

approach. 


E. TAPPERT SMOOTHING FUNCTION APPROACH 


In 1991, Tappert [13] suggested an alternative to the effective index term in the 
smoothing function as 


~2 2 
n=n 


a d^H[g) 
dz^ 


a = 


^-[p^lpbY 

1 


(44) 


This method can utilize the existing MMPE code by simply replacing the effective index 
of refraction term with Equation (44). We shall utilize the same definition for the mixing 
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length, L = —, although Tappert believed this should improve the model as both Az and 

K 

L get smaller. We shall not investigate the effect of reducing the mixing length with this 
approach further in this thesis. 
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III. RESULTS 


A. MMPE USER INTERFACE 

The MMPE input parameters give the user the ability to manually define the grid 
size. This is accomplished by changing the “pefiles.inp” input file in the local directory. 
The depth grid size is defined indirectly by changing nz in “pefiles.inp” to the desired 
length of the EFT size to be used in the calculation. This value determines the number of 
depth samples at each range step according to Az = '^Zj^jnz, where z^^^ is the maximum 
computational depth. The factor 2 accounts for the image solution method employed in 
SSF algorithms. Similarly, the range step size can be manually defined by changing rfact. 
This defines the range step by multiplying Az by rfact. Therefore, an rfact = 1 makes a 
uniform grid in which the depth grid size and the range step size are equal. An rfact >1 
makes the range step size larger than the depth grid size (and vice versa). If Az and rfact 
are set to 0.0, the model will automatically define Az = /I/IO and Ar = yi, where A, is the 
acoustic wavelength be computed. 

B. NO DENSITY CASE 

The original MMPE and both alternate options were used to model the 
transmission loss in an ocean environment with no density discontinuity at the bottom 
interface. In each case, the MMPE model was allowed to automatically determine the 
grid size based on the default values of Az and Ar, i.e., Az = /l/10 and Ar = /l. The 
solution for the alternate implementation methods were confirmed to be identical to the 
solution obtained in the original MMPE model, since there is no density discontinuity 
calculation in this example. All solutions closely match the benchmark solution, as 
observed in Figure 2. 

Since both alternate implementation options exactly match the original MMPE 
solution in this environment, only the original MMPE model was used to find the 
optimum Ar in the case without a density contrast. The Ar variability study was 
conducted by varying rfact while holding nz constant. 
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No density - accuracy 



Range (km) 

Figure 2. Original MMPE vs. Couple, no density discontinuity. 

a. Results form = 512 

Figures 3-6 show the results for a fixed value of nz = 512, corresponding to a Az 
of approximately 5 m, or Az~/l/3, with rfact values of 50, 20, 10, 5, 1, 0.2, 0.1, and 
0.05. These graphs show two notable features: that the model converges to a stable 
solution and that the stable solution tends to underestimate the transmission loss of 
the benchmark solution and lacks structure. Focusing on the transmission loss from 15- 
20 km, as depicted in Figure 7, allows us to better see the model’s convergence. For the 
nz = 512 case, the solution tends to converge at rfact values smaller than 1. 

Since an rfact = 5 overestimated transmission loss when compared with the 
benchmark solution, and rfact = 1 underestimated the transmission loss, there must be an 
optimum rfact between these values which provides the most accurate solution when 
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TL (dB re 1m) for Pressure 


compared to the benchmark. The optimum value for the nz = 5\2 case is determined to be 
2.5, as shown in Figure 8. This corresponds to a Ar of approximately 12 m, or /Sr ~ X. 


MMPE - No density - nz = 512 
” 1 “ 
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Figure 3. Original MMPE, no density, nz = 512, rfact = 50, 20. 


19 


























TL (dB re 1 m) for Pressure htI TL (dB re 1 m) for Pressure 


MMPE - No density - nz = 512 



tgure 4. 


Original MMPE, no density, nz = 512, rfact = 50, 20. 
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Figure 5. Original MMPE, no density, nz = 512, rfact = 1, 0.2. 
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Figure 6. Original MMPE, no density, nz = 512, rfact = 0.1, 0.05. 


MMPE - No density - nz = 512 



Figure 7. Solution convergence for nz = 512. 
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MMPE - No density - nz = 512 



Figure 8. Closest approximation to benchmark solution, nz = 512. 
b. Results for nz = 1024 

The process was repeated as Azwas decreased with similar results. Here, the 
depth mesh is approximately 10 m, or Az~ Aj6. The solution again converges at a value 
that slightly underestimates transmission loss as observed in Figures 9-13. In this case, 
the optimum rfact is determined to be 5 as shown in Figure 14. This again corresponds to 
a Ar of approximately 12 m, or Ar ~ A,. 
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Figure 9. Original MMPE, no density, nz = 1024, rfact = 50, 20. 
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Figure 10. Original MMPE, no density, nz = 1024, rfact = 10, 5. 
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Figure 11. Original MMPE, no density, nz = 1024, rfact = 1, 0.2. 
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Figure 12. Original MMPE, no density, nz = 1024, rfact = 0.1, 0.05. 
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MMPE - No density - nz = 1024 



Figure 13. Solution convergence for nz = 1024. 


MMPE - No density - nz = 1024 



Figure 14. Closest approximation to benchmark solution, nz = 1024. 
c. Results fornz>1024 


The model was also evaluated at nz values of 2048, 4096, 8192, and 16,384. The 
results for smaller Azwere similar to those obtained for nz = 512 and 1024. In all cases, 
the solution converges at a value slightly lower than the benchmark solution and the 
optimum solution occurs when Aris approximately 12 m, or Ar~/l. See the Appendix 
for full results. 
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c. 


MODEL RESULTS IN PRESENCE OF DENSITY DISCONTINUITY 


All three implementation options for the MMPE model were employed in an 
environment with a density contrast at the bottom interface. The range step size was held 
constant at 12 m while the depth grid size was reduced for each implementation method. 
Range step correction factors slightly higher and lower were also tested against a varying 
depth grid size. 

I. Original MMPE 

With a fixed Ar=12 m, there are negligible differences in the transmission loss 
model as observed in Figure 15. In this case, an nz value of 512 is not only the most 
efficient calculation, but also provides a transmission loss solution that is not 
significantly different than those obtained from a finer depth grid. As noted, for 
nz = 512, the depth mesh is approximately Az ~ 5 m ~ /1/3. 


Original MMPE - with density - deltar = 12 m 



Figure 15. Original MMPF, Ar= 12 m, various nz. 


2. MMPE with Tapper! Smoothing Function 

Figure 16 shows similar results when the Tappert smoothing function is 
implemented into the MMPF model. When compared to the original MMPF, the 
differences in the model are even less significant as Az is decreased. 
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MMPE with correction - with density - deltar = 12 m 



Figure 16. MMPE with Tappert smoothing function, Ar= 12 m, various nz. 

3. Yevick and Thomson Hybrid Implementation 

Yevick and Thomson’s SSF/FD hybrid implementation provides for some 
interesting results. At the largest Az (nz = 512), the model is a poor comparison to the 
benchmark solution at range. The solution underestimates the benchmark solution at all 
values. Furthermore the shape of the graph only vaguely represents the benchmark 
solution. Decreasing the depth grid size by a factor of 2 results in an even worse solution 
as seen in Figure 17. However, decreasing Az beyond this point improved the solution 
dramatically. At nz = 2048, the solution begins approaching the benchmark solution, both 
in overall shape and magnitude of transmission loss. The hybrid model most closely 
approximates the benchmark solution at nz = 8192, corresponding to Az « 0.3 m - 2/50 . 
Decreasing Az further results in smaller changes to the solution as observed in Figure 18. 
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SSF/FD hybrid - with density - deltar = 12 m 



Figure 17. Hybrid implementation, Ar= 12 m, various nz (coarse Az). 


SSF/FD hybrid - with density - deltar = 12 m 



Figure 18. Hybrid implementation, Ar= 12 m, various nz (fine Az). 


4. Overall Comparison 

Solutions for the original MMPE, MMPE with Tappert smoothing function, and 
Yevick and Thomson’s hybrid implementation are compared in Eigure 19. An nz value of 
512(Az~/l/3)is used for both the original version of the MMPE model and the MMPE 
model with the Tappert alternative method, since that value was efficient and smaller Az 
values did not significantly affect the solution accuracy. An nz value of 8192 (Az ~ /1/50) 
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is used for the hybrid implementation option since that value provided for the most 
accurate solution for the hybrid implementation technique. The figure shows there is 
marginal difference between the original MMPE and MMPE with the Tappert alternative 
function. The hybrid implementation option provides for a slightly closer approximation 
to the benchmark solution in amplitude, however is a less efficient solution. In all cases, 
though, the phase error exists at this long range. Thus, it does not appear that this phase 
error is caused by inaccuracies in the computational treatment of the smoothing function, 
or the first order correction in the hybrid method defined in Equation (22). 


Implementation Comparison - with density - deltar = 12 m 



Eigure 19. Implementation method comparison. 


5. Effect of Reference Sound Speed 

To this point, the reference sound speed for all transmission loss solutions was 
held constant at 1500 m/s. However, varying the reference sound speed has a significant 
effect on the phase of the solution. In general, the magnitudes of the solutions are 
unchanged, but the phase error present can be minimized. Eigures 20-23 show the results 
when the reference sound speed is varied. The previously determined optimum depth 
mesh and range step is used for each implementation option. The best solution appears to 
occur at a reference sound speed of 1480 m/s. Eigure 24 shows the results of all three 
implementation options at the optimum reference sound speed. 
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Original MM PE - with density - deltar = 12 m - various c 



Range (km) 


Figure 20. Original MMPE, Ar = 12 m, nz = 512, various Cq. 


MM PE with correction - with density - deltar = 12 m - various c 



Figure 21. MMPE with Tappert smoothing, Ar = 12 m, nz = 512, various Cq. 
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SSF/FD hybrid - with density - deltar = 12 m - various c 
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gure 22. 


Hybrid implementation, Ar = 12 m, nz = 8192, various Cq. 


Implementation comparison - with density - deltar = 12 m - c0=1480 



Figure 23. Implementation comparison, Ar 


12 m, 15-20 km. 
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Implementation comparison - with density- deltar = 12 m - c0=1480 
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Figure 24. Implementation comparison, Ar = 12 m, 0-20 km. 


The model’s inaccuracies at range appear primarily due to the phase errors. 
Alternate implementation options did not significantly correct the phase error. Proper 
selection of reference sound speed seems to correct this error. This was previously noted 
in a deep water benchmark problem defined in the 1993 PE Workshop [12], and 
subsequent work by Tappert [13] using his “co-insensitive” version seemed to correct the 
phase errors. Yevick and Thomson [10] were also able to correct the phase error found in 
the deep water PE Workshop problem by implementing their higher-order correction 
term, defined in Equation (38), in their hybrid method. It is quite possible that 
implementation of that higher-order correction term would also address the phase error 
found in the shallow water problem defined in this thesis. Although not as efficient as a 
basic SSE method, the proposed hybrid method with higher-order correction might be 
expected to improve the long-range phase accuracy of all solutions. 
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IV. SUMMARY 


The parabolic equation method for transmission loss modeling is an important 
topic in the acoustics scientific community. The SSF technique has been used for years 
due to its computational efficiency. Presently, SSF techniques provide accurate solutions 
in deep ocean environments where bottom interactions are negligible. However, as the 
focus shifts to shallow water acoustic propagation, it becomes necessary to find a way to 
use the SSF algorithm to model transmission loss. This paper compares two different 
implementation techniques when compared to a standard SSF model. 

The basic SSF algorithm employing the Thomson-Chapman WAPE 
approximation employed in the MMPE model was found to be very accurate when there 
is no density discontinuity present. As the depth mesh grid is reduced in size, the MMPE 
model converges to a stable solution. When compared to the benchmark solution, the 
convergence in Ar slightly underestimates transmission loss at all ranges. This is 
presumably due to numerical noise when the range step becomes too small. By 
systematically reducing Ar it is shown that an optimal comparison exists at a Ar ~ T . 
This optimal value for Ar remains constant even as Az is decreased and a density 
discontinuity is introduced. 

In the presence of a density discontinuity, both the original MMPE, and MMPE 
with the alternate Tappert approach to treating the effective index of refraction function, 
provide a fair model of the magnitude of transmission loss at ranges greater than 15 km. 
The largest tested Az of approximately /1/4 provided acceptable results. Decreasing Az 
further does not significantly change the solutions’ accuracy. Since a smaller Az results 
in increased computational time, it is recommended that a Az of approximately /1/4 be 
used for both options. 

The implementation of Yevick and Thomson’s hybrid technique was able to 
provide a slightly more accurate comparison to the benchmark solution amplitude. 
However, it is necessary to make Az very small in order to obtain such comparable 


33 



results. This comes at the cost of increased computational complexity, and should be 
balanced against the need to obtain the most accurate result. 

All models were found to suffer from a phase error, which accumulates in range 
when a reference sound speed of 1500 m/s is used. This phase error can be minimized by 
lowering cq to 1480 m/s. This is reminiscent of previous phase errors found with the 
Thomson-Chapman WAPE approximation. It is quite possible that these phase errors 
may be significantly reduced by employing the higher-order correction of Yevick and 
Thomson, or the Co-insensitive version of the SSF algorithm defined by Tappert. Future 
work could investigate the relative accuracy of each of these methods. 
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APPENDIX. MMPE CONVERGENCE GRAPHS 
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Figure 25. Original MMPE, no density, nz = 2048, rfact = 50, 20. 
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Figure 26. Original MMPE, no density, nz = 2048, rfact = 10, 5. 
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Figure 27. Original MMPE, no density, nz = 2048, rfact = 1, 0.2. 


MMPE-No demity-nz - 2048 



Figure 28. Original MMPE, no density, nz = 2048, rfact = 0.1, 0.05. 
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MMPE - No density - nz = 2048 



Range (km) 


Figure 29. Solution convergence for nz = 2048. 


MMPE - No density - nz = 2048 



Figure 30. Closest approximation to benchmark solution, nz = 2048. 
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Figure 31. Original MMPE, no density, nz = 4096, rfact = 100, 50. 
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Figure 32. Original MMPE, no density, nz = 4096, rfact = 20, 10. 
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Figure 33. Original MMPE, no density, nz = 4096, rfact = 5, 1. 


MMPE - No density - nz = 4096 



Figure 34. Solution convergence for nz = 4096. 
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MMPE - No density - nz = 4096 



Figure 35. Closest approximation to benchmark solution, nz = 4096. 
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Figure 36. Original MMPE, no density, nz = 8192, rfact = 200, 100. 
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Figure 37. Original MMPE, no density, nz = 8192, rfact = 50, 20. 


MMPE-No denoKy-nz- 8192 



Figure 38. Original MMPE, no density, nz = 8192, rfact = 10, 5. 
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MMPE - No density - nz = 8192 
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Figure 39. Solution convergence for nz = 8192. 


MMPE - No density - nz = 8192 



Figure 40. Closest approximation to benchmark solution, nz = 8192. 
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Figure 41. Original MMPE, no density, nz = 16,384, rfact = 800, 400. 
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Figure 42. Original MMPE, no density, nz = 16,384, rfact = 200, 100. 
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Figure 43. Original MMPE, no density, nz = 16,384, rfact = 50, 20. 
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Figure 44. Solution convergence for nz = 16,384. 
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MMPE - No density - nz = 16384 



Figure 45. Closest approximation to benchmark solution, nz = 16,384. 
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